{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Likelihood Analysis with Python\n",
    "\n",
    "The python likelihood tools are a very powerful set of analysis tools that expand upon the command line tools provided with the Fermitools package.\n",
    "\n",
    "Not only can you perform all of the same likelihood analysis with the python tools that you can with the standard command line tools, but you can directly access all of the model parameters.\n",
    "\n",
    "You can more easily script a standard analysis like light curve generation. There are also a few things built into the python tools that are not available from the command line like the calculation of upper limits.\n",
    "\n",
    "There are many [user contributed packages](https://fermi.gsfc.nasa.gov/ssc/data/analysis/user/) built upon the python backbone of the Fermitools and we are going to highlight and use a few of them in this tutorial like `likeSED`, `make4FGLxml`, and the `LATAnalysisScripts`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This sample analysis is based on the PG 1553+113 analysis performed by the LAT team and described in [Abdo, A. A. et al. 2010, ApJ, 708, 1310](http://adsabs.harvard.edu/abs/2010ApJ...708.1310A).\n",
    "\n",
    "At certain points we will refer to this article as well as the [Cicerone](https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/).\n",
    "\n",
    "After you complete this tutorial, you should be able to reproduce all of the data analysis performed in this publication including generating a spectrum (individual bins and a butterfly plot) and produce a light curve with the python tools.\n",
    "\n",
    "This tutorial assumes you have the most recent [Fermitools](https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/) installed.\n",
    "\n",
    "We will also make significant use of python, so you might want to familiarize yourself with python (there's a beginner's guide at http://wiki.python.org/moin/BeginnersGuide). This tutorial also assumes that you've gone through the non-python based unbinned likelihood thread. This tutorial should take approximately 8 hours to complete (depending on your computer's speed) if you do everything, but there are some steps you can skip along the way which shave about 4 hours off of that."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get the data\n",
    "\n",
    "For this thread the original data were extracted from the [LAT data server](https://fermi.gsfc.nasa.gov/cgi-bin/ssc/LAT/LATDataQuery.cgi) with the following selections (these selections are similar to those in the paper):\n",
    "\n",
    "```\n",
    "Search Center (RA,Dec) = (238.929,11.1901)\n",
    "Radius = 20 degrees\n",
    "Start Time (MET) = 239557417 seconds (2008-08-04T15:43:37)\n",
    "Stop Time (MET) = 256970880 seconds (2009-02-22T04:48:00)\n",
    "Minimum Energy = 100 MeV\n",
    "Maximum Energy = 300000 MeV\n",
    "```\n",
    "\n",
    "For more information on how to download LAT data please see the [Extract LAT Data](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/extract_latdata.html) tutorial.\n",
    "\n",
    "These are the event files. Run the code cell below to retrieve them:\n",
    "```\n",
    "L1504241622054B65347F25_PH00.fits\n",
    "L1504241622054B65347F25_PH01.fits\n",
    "L1504241622054B65347F25_SC00.fits\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/pyLikelihood/L1504241622054B65347F25_PH00.fits\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/pyLikelihood/L1504241622054B65347F25_PH01.fits\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/pyLikelihood/L1504241622054B65347F25_SC00.fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mkdir data\n",
    "!mv *SC00.fits PG1553_SC.fits\n",
    "!mv *.fits ./data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You'll first need to make a file list with the names of your input event files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls ./data/*PH*.fits > ./data/PG1553_events.list\n",
    "!cat ./data/PG1553_events.list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following analysis, we've assumed that you've named your list of data files `PG1553_events.list` and the spacecraft file `PG1553_SC.fits`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Perform Event Selections\n",
    "\n",
    "You could follow the unbinned likelihood tutorial to perform your event selections using **gtlike*, **gtmktime**, etc. directly from the command line, and then use pylikelihood later.\n",
    "\n",
    "But we're going to go ahead and use python. The `gt_apps` module provides methods to call these tools from within python. This'll get us used to using python.\n",
    "\n",
    "So, let's jump into python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gt_apps as my_apps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, you can see what objects are part of the gt_apps module by executing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(my_apps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The python object for **gtselect** is called `filter`, and we first need to set all of its options.\n",
    "\n",
    "This is very similar to calling **gtselect** from the command line and inputting all of the variables interactively. It might not seem that convenient to do it this way, but it's really nice once you start building up scripts (see [Building a Light Curve](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/python_tutorial.html#lightcurve)) and reading back options and such.\n",
    "\n",
    "For example, towards the end of this thread, we'll want to generate a light curve and we'll have to run the likelihood analysis for each datapoint. It'll be much easier to do all of this within python and change the `tmin` and `tmax` in an iterative fashion.\n",
    "\n",
    "Note that these python objects are just wrappers for the standalone tools, so if you want any information on their options, see the corresponding [documentation](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/references.html) for the standalone tool.\n",
    "\n",
    "Let's set **gtselect**'s (or `filter`'s) options:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.filter['evclass'] = 128\n",
    "my_apps.filter['evtype'] = 3\n",
    "my_apps.filter['ra'] = 238.929\n",
    "my_apps.filter['dec'] = 11.1901\n",
    "my_apps.filter['rad'] = 10\n",
    "my_apps.filter['emin'] = 100\n",
    "my_apps.filter['emax'] = 300000\n",
    "my_apps.filter['zmax'] = 90\n",
    "my_apps.filter['tmin'] = 239557417\n",
    "my_apps.filter['tmax'] = 256970880\n",
    "my_apps.filter['infile'] = '@./data/PG1553_events.list'\n",
    "my_apps.filter['outfile'] = './data/PG1553_filtered.fits'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now we run **gtselect**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.filter.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that you can see exactly what gtselect will do if you run it by typing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.filter.command()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You have access to any of the inputs by directly accessing the `filter['OPTIONS']` options.\n",
    "\n",
    "Next, you need to run **gtmktime**. This is accessed within python via the `maketime` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.maketime['scfile'] = './data/PG1553_SC.fits'\n",
    "my_apps.maketime['filter'] = '(DATA_QUAL>0)&&(LAT_CONFIG==1)'\n",
    "my_apps.maketime['roicut'] = 'no'\n",
    "my_apps.maketime['evfile'] = './data/PG1553_filtered.fits'\n",
    "my_apps.maketime['outfile'] = './data/PG1553_filtered_gti.fits'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.maketime.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're using the most conservative and most commonly used cuts described in deatil in the [Cicerone](http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Likelihood/Exposure.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Livetime Cubes and Exposure Maps\n",
    "\n",
    "At this point, you could make a counts map of the events we just selected using [gtbin](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtbin.txt) (it's called `evtbin` within python) and I won't discourage you, but we're going to go ahead and create a livetime cube and exposure map.\n",
    "\n",
    "This might take a few minutes to complete, so if you want to create a counts map and have a look at it, get these processes going and open another terminal to work on your counts map (see the [likelihood tutorial](http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/likelihood_tutorial.html#makeCountsMaps) for an example of running **gtbin** to produce a counts map)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Livetime Cube"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This step will take approximately 15 - 30 minutes to complete so if you want to just download the `PG1553_ltCube.fits` from us you can skip this step by running the code cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/pyLikelihood/PG1553_ltCube.fits\n",
    "!mv *.fits ./data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If instead you wish to do this yourself, **gtltcube** is `expCube` and you set its options like so:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.expCube['evfile'] = './data/PG1553_filtered_gti.fits'\n",
    "my_apps.expCube['scfile'] = './data/PG1553_SC.fits'\n",
    "my_apps.expCube['outfile'] = './data/PG1553_ltCube.fits'\n",
    "my_apps.expCube['zmax'] = 90\n",
    "my_apps.expCube['dcostheta'] = 0.025\n",
    "my_apps.expCube['binsz'] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.expCube.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While you're waiting, you might have noticed that not all of the command line Fermitools have an equivalent object in `gt_apps`.\n",
    "\n",
    "This is easy to fix. Say you want to use **gtltcubesun** from within python. Just make it a GtApp:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from GtApp import GtApp\n",
    "expCubeSun = GtApp('gtltcubesun','Likelihood')\n",
    "expCubeSun.command()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exposure Map\n",
    "\n",
    "Similar to Livetime Cube above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.expMap['evfile'] = './data/PG1553_filtered_gti.fits'\n",
    "my_apps.expMap['scfile'] = './data/PG1553_SC.fits'\n",
    "my_apps.expMap['expcube'] = './data/PG1553_ltCube.fits'\n",
    "my_apps.expMap['outfile'] = './data/PG1553_expMap.fits'\n",
    "my_apps.expMap['irfs'] = 'CALDB'\n",
    "my_apps.expMap['srcrad'] = 20\n",
    "my_apps.expMap['nlong'] = 120\n",
    "my_apps.expMap['nlat'] = 120\n",
    "my_apps.expMap['nenergies'] = 37"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.expMap.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate XML Model File\n",
    "\n",
    "We need to create an XML file with all of the sources of interest within the Region of Interest (ROI) of PG 1553+113 so we can correctly model the background. For more information on the format of the model file and how to create one, see the [likelihood analysis](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/likelihood_tutorial.html#createSourceModel) tutorial.\n",
    "\n",
    "For this thread, we use a simplified model with just two sources: `PG1553compact_model.xml`.\n",
    "\n",
    "To fit all `4FGL` sources, you can use the user-contributed tool `make4FGLxml.py` to create a model file based on the LAT 8-year LAT catalog. You'll need to download the FITS version of this file at http://fermi.gsfc.nasa.gov/ssc/data/access/lat/8yr_catalog/ and get the `make4FGLxml.py` tool from the [user-contributed software page](https://fermi.gsfc.nasa.gov/ssc/data/analysis/user/) and put them both in your working directory. Also make sure you have the most recent galactic diffuse and isotropic model files which can be found at http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html. \n",
    "\n",
    "The catalog and background models are also packaged with your installation of the Fermitools, which can be found at: `$FERMI_DIR//refdata/fermi/galdiffuse/`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/user/make4FGLxml.py\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/access/lat/8yr_catalog/gll_psc_v19.fit\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/pyLikelihood/PG1553compact_model.xml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mv *.fit *.xml ./data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once you have all of the files, you can generate your model file in python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from make4FGLxml import *\n",
    "mymodel = srcList('./data/gll_psc_v19.fit','./data/PG1553_filtered_gti.fits','./data/PG1553_model.xml')\n",
    "mymodel.makeModel('./data/gll_iem_v07.fits', 'gll_iem_v07', './data/iso_P8R3_SOURCE_V2_v1.txt', 'iso_P8R3_SOURCE_V2_v1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more information on the `make4FGLxml.py` module, see the [usage notes](https://fermi.gsfc.nasa.gov/ssc/data/analysis/user/readme_make4FGLxml.txt).\n",
    "\n",
    "In the paper, the LAT team only included two sources: one from the `0FGL` catalog and another non-catalog source. This is because the later LAT catalogs had not been released at the time. For a more updated model, one should use the latest 4FGL information.\n",
    "\n",
    "Back to looking at our compact XML model file: notice that PG 1553+113 is listed in the model file as 4FGL J1555.7+1111 with all of the parameters filled in for us.\n",
    "\n",
    "It's actually offset from the center of our ROI by 0.003 degrees.\n",
    "\n",
    "How nice!\n",
    "\n",
    "Also notice that we changed the model for 4FGL J1555.7+1111 from a Log Parabola to a simple power-law for the purposes of this analysis thread."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute the diffuse source responses\n",
    "\n",
    "The diffuse source responses tell the likelihood fitter what the expected contribution would be for each diffuse source, given the livetime associated with each event.\n",
    "\n",
    "The source model XML file must contain all of the diffuse sources to be fit. The [gtdiffrsp](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtdiffrsp.txt) tool will add one column to the event data file for each diffuse source.\n",
    "\n",
    "The diffuse response depends on the instrument response function (IRF), which must be in agreement with the selection of events, i.e. the event class and event type we are using in our analysis. Since we are using SOURCE class, `CALDB` should use the `P8R3_SOURCE_V2` IRF for this tool.\n",
    "\n",
    "If the diffuse responses are not precomputed using [gtdiffrsp](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtdiffrsp.txt), then the [gtlike](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtlike.txt) tool will compute them at runtime (during the next step).\n",
    "\n",
    "However, as this step is very computationally intensive (often taking ~hours to complete), and it is very likely you will need to run [gtlike](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtlike.txt) more than once, it is probably wise to precompute these quantities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.diffResps['evfile'] = './data/PG1553_filtered_gti.fits'\n",
    "my_apps.diffResps['scfile'] = './data/PG1553_SC.fits'\n",
    "my_apps.diffResps['srcmdl'] = './data/PG1553compact_model.xml'\n",
    "my_apps.diffResps['irfs'] = 'CALDB'\n",
    "my_apps.diffResps.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run the Likelihood Analysis\n",
    "\n",
    "It's time to actually run the likelihood analysis now.\n",
    "\n",
    "First, you need to import the pyLikelihood module and then the UnbinnedAnalysis functions (there's also a binned analysis module that you can import to do [binned likelihood analysis](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/binned_likelihood_tutorial.html) which behaves almost exactly the same as the unbinned analysis module).\n",
    "\n",
    "For more details on the pyLikelihood module, check out the [pyLikelihood Usage Notes](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/python_usage_notes.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyLikelihood\n",
    "from UnbinnedAnalysis import *\n",
    "\n",
    "obs = UnbinnedObs('./data/PG1553_filtered_gti.fits','./data/PG1553_SC.fits',expMap='./data/PG1553_expMap.fits', expCube='./data/PG1553_ltCube.fits',irfs='CALDB')\n",
    "like = UnbinnedAnalysis(obs,'./data/PG1553compact_model.xml',optimizer='NewMinuit')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By now, you'll have two objects: `obs`, an UnbinnedObs object, and `like`, an UnbinnedAnalysis object.\n",
    "\n",
    "You can view these objects attributes and set them from the command line in various ways. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(like)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or you can get directly at the objects attributes and methods by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir(like)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or get even more details by executing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(like)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are a lot of attributes, and here you start to see the power of using pyLikelihood since you'll be able to (once the fit is done) access any of these attributes directly within python and use them in your own scripts.\n",
    "\n",
    "For example: you can see that the like object has a `tol` attribute, which we can read back to see what it is and, then set it to what we want it to be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "like.tol\n",
    "like.tolType\n",
    "like.tol = 0.0001"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we're ready to do the actual fit.\n",
    "\n",
    "This next step can take anywhere from 10 minutes to a few hours to complete depending on your system.\n",
    "\n",
    "We're doing something a bit fancy here. We're getting the minimizating object (and calling it `likeobj`) from the `logLike` object so that we can access it later.\n",
    "\n",
    "We pass this object to the fit routine so that it knows which fitting object to use. We're also telling the code to calculate the covariance matrix so we can get at the errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "likeobj = pyLike.NewMinuit(like.logLike)\n",
    "like.fit(verbosity=0,covar=True,optObject=likeobj)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number that is printed out here is the -log(Likelihood) of the total fit to the data.\n",
    "\n",
    "You can print the results of the fit by accessing `like.model`. You can now plot the results of the fit by executing the `plot` command. The results are shown below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "like.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output of the plot function of the like1 UnbinnedAnalysis object shows:\n",
    "\n",
    "* Left: the contribution of each of the objects in the model to the total model, and plots the data points on top.\n",
    "\n",
    "* Right: the residuals of the likelihood fit to the data. Notice that the fit is poor in the second to last bin.\n",
    "\n",
    "Now, check if `NewMinuit` converged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print likeobj.getRetCode()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you get anything other than '0' here, then NewMinuit didn't converge. There are several reasons that we might not get convergence:\n",
    "\n",
    "* The culprit is usually a parameter (or parameters) in the model that reach the limits set in the XML file. If this happens, the minimizer cannot reach the formal minimum and hence cannot calculate the curvature.\n",
    "\n",
    "\n",
    "* Often the problem is with spectral shape parameters (PL index, etc.), so simply freezing the shape of all spectral parameters to their values from the full time period (and certainly for weaker background sources) when fitting a shorter time period may solve the problem. Remember that the 4FGL catalog used a full 8 years of data and we're using a much shorter time period here.\n",
    "\n",
    "\n",
    "* Weak background sources are more likely to cause problems, so you could consider just freezing them completely (or removing them from the model).\n",
    "\n",
    "    For example: a background source from the catalog that is detected at TS~=25 in 2 years could cause convergence problems in a 1-month light curve, where it will often not be detectable at all.\n",
    "\n",
    "\n",
    "* If there are no parameters at their limits, then increasing the overall convergence tolerance may help - try using a value of 1E-8 for the absolute tolerance.\n",
    "\n",
    "\n",
    "* If that doesn't help, then try to systematically simplify the model. Progressively freeze all sources, starting with those at the edge of the ROI in and moving in until you get a model simple enough for the minimizer to work reliably.\n",
    "\n",
    "    For example: if you are using a 10 degree ROI, you could start by freezing all background sources further than 7 degrees from the source of interest, and move to 5 degrees if that doesn't solve the problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our case, we have convergence, and we can move forward.\n",
    "\n",
    "Note, however, that this doesn't mean you have the 'right' answer. It just means that you have the answer assuming the model you put in. This is a subtle feature of the likelihood method.\n",
    "\n",
    "Let's take a look at the residuals and overall model using `matplotlib`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to import these two modules for plotting and some other stuff. They are included in the Fermitools."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = (like.energies[:-1] + like.energies[1:])/2.\n",
    "# The 'energies' array are the endpoints so we take\n",
    "# the midpoint of the bins.\n",
    "\n",
    "plt.figure(figsize=(5,5))\n",
    "plt.ylim((0.4,1e4))\n",
    "plt.xlim((200,300000))\n",
    "sum_model = np.zeros_like(like._srcCnts(like.sourceNames()[0]))\n",
    "\n",
    "for sourceName in like.sourceNames():\n",
    "    sum_model = sum_model + like._srcCnts(sourceName)\n",
    "    plt.loglog(E,like._srcCnts(sourceName),label=sourceName[1:])\n",
    "\n",
    "plt.loglog(E,sum_model,label='Total Model')\n",
    "plt.errorbar(E,like._Nobs(),yerr=np.sqrt(like._Nobs()), fmt='o',label='Counts')\n",
    "ax = plt.gca()\n",
    "box = ax.get_position()\n",
    "ax.set_position([box.x0, box.y0, box.width * 0.5, box.height])\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resid = (like._Nobs() - sum_model)/sum_model\n",
    "resid_err = (np.sqrt(like._Nobs())/sum_model)\n",
    "plt.figure(figsize=(9,9))\n",
    "plt.xscale('log')\n",
    "plt.errorbar(E,resid,yerr=resid_err,fmt='o')\n",
    "plt.axhline(0.0,ls=':')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are just some examples of how to access the underlying data from the python tools.\n",
    "\n",
    "As you can see, most of the information is accessable if you dig around. Note that the residuals aren't that great in several bins especially at high energies.\n",
    "\n",
    "This could be a missing source or that we're not modelling a bright source very well. If you look up at the model plot, the high energy tail of PG1553 might be biasing the fit at high energies and we might get a better fit with a log-parabola or broken power-law. That's for another day, though.\n",
    "\n",
    "Let's check out the final parameters of the fit for PG1553."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "like.model['4FGL J1555.7+1111']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "like.flux('4FGL J1555.7+1111',emin=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "like.fluxError('4FGL J1555.7+1111',emin=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also get the TS value for a specific source:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "like.Ts('4FGL J1555.7+1111')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can calculate how many standard deviations (σ) this corresponds to.\n",
    "\n",
    "Remember that the [TS value is ∼ $σ^2$](https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Likelihood/Likelihood_overview.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sqrt(like.Ts('4FGL J1555.7+1111'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "like.logLike.writeXml('PG1553_fit.xml')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mv *.xml ./data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create some TS Maps\n",
    "\n",
    "If you want to check on your fit, you should probably make a test statistic (TS) map of your region at this point.\n",
    "\n",
    "This will allow you to check and make sure that you didn't miss anything in the ROI (like a source not in the catalog you used).\n",
    "\n",
    "A TS map is created by moving a putative point source through a grid of locations on the sky and maximizing the log(likelihood) at each grid point.\n",
    "\n",
    "We're going to make two TS Maps. If you have access to multiple CPUs, you might want to start one of the maps running in one terminal and the other in another one since this will take **several hours** to complete.\n",
    "\n",
    "Our main goal is to create one TS Map with PG1553 included in the fit and another without. The first one will allow us to make sure we didn't miss any sources in our ROI and the other will allow us to see the source.\n",
    "\n",
    "All of the other sources in our model file will be included in the fit and shouldn't show up in the final map."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you don't want to create these maps, go ahead and skip the following steps and jump to the section on making a butterfly plot. The creation of these maps is optional and is only for double checking your work.\n",
    "\n",
    "You should still modify your model fit file as detailed in the following:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copy `PG1553_fit.xml` as `PG1553_fit_TSMap.xml`. Modify `PG1553_fit_TSMap.xml` to make everything fixed by changing all of the `free=\"1\"` statements to `free=\"0\"`. Modify the lines that look like:\n",
    "\n",
    "```xml\n",
    "<parameter error=\"1.670710816\" free=\"1\" max=\"10000\" min=\"0.0001\" name=\"Prefactor\" scale=\"1e-13\" value=\"6.248416874\" />\n",
    "<parameter error=\"0.208035277\" free=\"1\" max=\"5\" min=\"0\" name=\"Index\" scale=\"-1\" value=\"2.174167564\" />\n",
    "```\n",
    "\n",
    "to look like:\n",
    "\n",
    "```xml\n",
    "<parameter error=\"1.670710816\" free=\"0\" max=\"10000\" min=\"0.0001\" name=\"Prefactor\" scale=\"1e-13\" value=\"6.248416874\" />\n",
    "<parameter error=\"0.208035277\" free=\"0\" max=\"5\" min=\"0\" name=\"Index\" scale=\"-1\" value=\"2.174167564\" />\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cp ./data/PG1553_fit.xml ./data/PG1553_fit_TSMap.xml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"./data/PG1553_fit_TSMap.xml\") as file:\n",
    "    file.read().replace(\"free=0\", \"free=1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then copy it and call it `PG1553_fit_noPG1553.xml` and comment out (or delete) the PG1553 source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cp ./data/PG1553_fit_TSMap.xml ./data/PG1553_fit_noPG1553.xml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " If you want to generate the TS maps, in one window do the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.TsMap['statistic'] = \"UNBINNED\"\n",
    "my_apps.TsMap['scfile'] = \"./data/PG1553_SC.fits\"\n",
    "my_apps.TsMap['evfile'] = \"./data/PG1553_filtered_gti.fits\"\n",
    "my_apps.TsMap['expmap'] = \"./data/PG1553_expMap.fits\"\n",
    "my_apps.TsMap['expcube'] = \"./data/PG1553_ltCube.fits\"\n",
    "my_apps.TsMap['srcmdl'] = \"./data/PG1553_fit_TSMap.xml\"\n",
    "my_apps.TsMap['irfs'] = \"CALDB\"\n",
    "my_apps.TsMap['optimizer'] = \"NEWMINUIT\"\n",
    "my_apps.TsMap['outfile'] = \"./data/PG1553_TSmap_resid.fits\"\n",
    "my_apps.TsMap['nxpix'] = 25\n",
    "my_apps.TsMap['nypix'] = 25\n",
    "my_apps.TsMap['binsz'] = 0.5\n",
    "my_apps.TsMap['coordsys'] = \"CEL\"\n",
    "my_apps.TsMap['xref'] = 238.929\n",
    "my_apps.TsMap['yref'] = 11.1901\n",
    "my_apps.TsMap['proj'] = 'STG'\n",
    "my_apps.TsMap.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In another window, do the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.TsMap['statistic'] = \"UNBINNED\"\n",
    "my_apps.TsMap['scfile'] = \"./data/PG1553_SC.fits\"\n",
    "my_apps.TsMap['evfile'] = \"./data/PG1553_filtered_gti.fits\"\n",
    "my_apps.TsMap['expmap'] = \"./data/PG1553_expMap.fits\"\n",
    "my_apps.TsMap['expcube'] = \"./data/PG1553_ltCube.fits\"\n",
    "my_apps.TsMap['srcmdl'] = \"./data/PG1553_fit_noPG1553.xml\"\n",
    "my_apps.TsMap['irfs'] = \"CALDB\"\n",
    "my_apps.TsMap['optimizer'] = \"NEWMINUIT\"\n",
    "my_apps.TsMap['outfile'] = \"./data/PG1553_TSmap_noPG1553.fits\"\n",
    "my_apps.TsMap['nxpix'] = 25\n",
    "my_apps.TsMap['nypix'] = 25\n",
    "my_apps.TsMap['binsz'] = 0.5\n",
    "my_apps.TsMap['coordsys'] = \"CEL\"\n",
    "my_apps.TsMap['xref'] = 238.929\n",
    "my_apps.TsMap['yref'] = 11.1901\n",
    "my_apps.TsMap['proj'] = 'STG'\n",
    "my_apps.TsMap.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This will take a while, so go get a cup of coffee.\n",
    "\n",
    "The final TS maps will be pretty rough since we selected to only do 25x25 0.5 degree bins, but it will allow us to check if anything is missing and go back and fix it.\n",
    "\n",
    "These rough plots aren't publication ready, but they'll do in a pinch to verify that we've got the model right."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the TS maps have been generated, open them up in whatever viewing program you wish.\n",
    "\n",
    "Here's how you could use `matplotlib` to do it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import astropy.io.fits as pyfits\n",
    "residHDU = pyfits.open('PG1553_TSmap_resid.fits')\n",
    "sourceHDU = pyfits.open('PG1553_TSmap_noPG1553.fits')\n",
    "fig = plt.figure(figsize=(16,8))\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(residHDU[0].data)\n",
    "plt.colorbar()\n",
    "plt.subplot(1,2,2)\n",
    "plt.imshow(sourceHDU[0].data)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that the `noPG1553` map (right) has a prominent source in the middle of the ROI.\n",
    "\n",
    "This is PG 1553+113, and if you execute the following, you can see that the TS value of the maxiumum pretty well with the TS value of PG 1553+113 that we arrived at above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.max(sourceHDU[0].data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you then look at the `residuals` plot (left), you can see that it's pretty flat in TS space.\n",
    "\n",
    "This indicates that we didn't miss any significant sources. If you did see something significant, you probably should go back and figure out what source you missed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Produce the Butterfly Plot\n",
    "\n",
    "Now we'll produce the butterfly plot.\n",
    "\n",
    "First, get the differential flux, the decorrelation (pivot) energy, and the index. You'll also need the errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N0 = like.model['4FGL J1555.7+1111'].funcs['Spectrum'].getParam('Prefactor').value()\n",
    "N0_err = like.model['4FGL J1555.7+1111'].funcs['Spectrum'].getParam('Prefactor').error()\n",
    "gamma = like.model['4FGL J1555.7+1111'].funcs['Spectrum'].getParam('Index').value()\n",
    "gamma_err = like.model['4FGL J1555.7+1111'].funcs['Spectrum'].getParam('Index').error()\n",
    "E0 = like.model['4FGL J1555.7+1111'].funcs['Spectrum'].getParam('Scale').value()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate the Butterfly Plot\n",
    "\n",
    "We need to calculate the differential flux at several energy points as well as the error on that differential flux. The differential flux is defined as:\n",
    "\n",
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/pyLikelihood/PG1553_thread8x.png'>\n",
    "\n",
    "and the error on that flux is defined as:\n",
    "\n",
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/pyLikelihood/PG1553_thread9x.png'>\n",
    "\n",
    "Where cov<sub>γγ</sub> = σγ<sup>2</sup>, where γ is the spectral index. To calculate cov<sub>γγ</sub>, use this code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "freeParValues = []\n",
    "\n",
    "for sourcename in like.sourceNames():\n",
    "    for element in like.freePars(sourcename):\n",
    "        freeParValues.append(element.getValue())\n",
    "\n",
    "g_index = freeParValues.index(like.freePars('4FGL J1555.7+1111')[1].getValue())\n",
    "cov_gg = like.covariance[g_index][g_index]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, let calculate F(E) and F(E) +/- Δ F(E) so that we can plot it.\n",
    "\n",
    "First, define a function for the flux and one for the error on the flux:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = lambda E,N0,E0,gamma: N0*(E/E0)**(-1*gamma)\n",
    "ferr = lambda E,F,N0,N0err,E0,cov_gg: F*np.sqrt(N0err**2/N0**2 + ((np.log(E/E0))**2)*cov_gg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, generate some energies to evaluate the functions above and evaluate them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = np.logspace(2,5,100)\n",
    "F = f(E,N0,E0,gamma)\n",
    "Ferr = ferr(E,F,N0,N0_err,E0,cov_gg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you can plot this. Here we're multiplying the flux by E<sup>2</sup> for clarity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,8))\n",
    "plt.xlabel('Energy [MeV]')\n",
    "plt.ylabel(r'E$^2$ dN/dE [MeV s$^{-1}$ cm$^{-2}$]')\n",
    "plt.loglog(E,E**2*(F+Ferr))\n",
    "plt.loglog(E,E**2*(F-Ferr))\n",
    "plt.plot(E,E**2*F)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The green line is the nominal fit, and the blue and orange lines are the 1 sigma contours."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate Spectral Points\n",
    "\n",
    "**Note**: you need a version of the Fermitools with ROOT installed to run through the steps in this section.\n",
    "\n",
    "To generate spectral points to plot on top of the butterfly that we just produced, you need to go back to the data selection part and use gtselect (filter in python) to divide up your data set in energy bins and run the likelihood fit on each of these individual bins.\n",
    "\n",
    "Luckily, there's a script that does this for us that we'll employ to get a final result. This script also generates the butterfly plot we just produced so you won't have to redo that again. It also does a lot of checks on the data to make sure that everything is going okay.\n",
    "\n",
    "If you've made it this far, you might be a little curious as to why we didn't jump right into using this tool, but now you're in a position to play with the python tools and make them do what you want them to.\n",
    "\n",
    "The script is also much more careful in handling units and saving everything to files than we have been in this interactive session."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So download the [likeSED.py](https://fermi.gsfc.nasa.gov/ssc/data/analysis/user/) user-contributed tool (it's in the SED_scripts package) and load it up into python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/user/SED_scripts_v13.1.tgz\n",
    "!mv *.tgz ./data\n",
    "!tar -xzvf ./data/SED_scripts_v13.1.tgz -C ./data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can find information on the usage of this tool on the [same page](https://fermi.gsfc.nasa.gov/ssc/data/analysis/user/) where you can download it.\n",
    "\n",
    "It can be used to generate an SED plot for both binned and unbinned analyses but we're only going to work on a binned analysis here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.insert(0, './data/SED_scripts_v13.1/')\n",
    "\n",
    "from likeSED import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Due to the SED scripts being slightly outdated, you may need to change a few things in your likeSED.py for it to import properly:\n",
    "\n",
    "* change `import pyfits` -> `import astropy.io.fits as pyfits`\n",
    "* change all instances of `pyfits.new_table` -> `pyfits.BinTableHDU.from_columns`, keeping all parameters intact."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you need to create a `likeInput` object that takes our unbinnedAnalysis object, the source name, and the number of bins we want as arguments.\n",
    "\n",
    "We're going to make 9 bins here like in the paper, and we're also going to make custom bins and bin centers. You can have the module chose bins and bin centers for you (via the `getECent` function), but we're going to do it this way so we can have the double-wide bin at the end.\n",
    "\n",
    "We're also going to use the `fit2` file (the result of our fit on the full dataset with the NewMinuit optimizer) for the fits of the individual bins, but we need to edit it first to make everything fixed except for the integral value of PG 1553+113.\n",
    "\n",
    "We're basically assuming that the spectral index is the same for all bins and that the contributions of the other sources within the ROI are not changing with energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inputs = likeInput(like,'4FGL J1555.7+1111',nbins=9)\n",
    "inputs.plotBins()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This last step will take a while (approximately 30 minutes) because it's creating an expmap and event file for each of the bins that we requested (look in your working directory for these files).\n",
    "\n",
    "Once it is done, we'll tell it to do the fit of the full energy band and make sure we request that it keep the covariance matrix.\n",
    "\n",
    "We will then create a likeSED object from our inputs, tell it where we want our centers to be, and then fit each of the bands. After this is done, we can plot it all with the Plot function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inputs.fullFit(CoVar=True)b\n",
    "sed = likeSED(inputs)\n",
    "sed.getECent()\n",
    "sed.fitBands()\n",
    "sed.Plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This produces the following three plots: a TS plot, a counts spectrum, and a butterfly plot (statistical errors only):\n",
    "\n",
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/pyLikelihood/likeSED_TSvEnergy.png'>\n",
    "\n",
    "The TS plot shows the TS value for each bin in the spectrum produced by likeSED. Notice that the last bin is double-wide.\n",
    "\n",
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/pyLikelihood/likeSED_cntSpec.png'>\n",
    "\n",
    "Counts spectrum produced by likeSED along with the spectral points.\n",
    "\n",
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/pyLikelihood/likeSED_EsqdNdE.png'>\n",
    "\n",
    "Butterfly plot along with the spectral points produced by likeSED. Note that this plot does not include the systematic errors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mmm, looking at that spectrum, there might be some curvature there like was hinted at in the model plot above. The butterfly model is the likelihood output and will only tell us the best fit power law since that is what we had in our model. The SED points are fits with a power-law in individual energy bands and can show us if the model we're assuming is a good assumption. Testing for curvature is for another day..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're going to leave python and look at the files that likeSED produced.\n",
    "\n",
    "Open a new terminal. In your working directory you will find a bunch of `4FGL_J1555.7+1111_9bins` files.\n",
    "\n",
    "These are just the event and exposure maps for each bin we selected. If you wanted to completely redo the likeSED analysis you should delete these.\n",
    "\n",
    "Some other files to note are `4FGL_J1555.7+1111_9bins_likeSEDout.fits` which contains all of the data needed to reproduce the plots seen when we ran the Plot function. There is also a root file (`4FGL_J1555.7+1111_9bins_histos.root`) with dumps of the plots and eps files of all of the plots.\n",
    "\n",
    "Now that we've generated a butterfly SED and individual SED data points, we can move on to working on a light curve."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Building a Light Curve\n",
    "\n",
    "TBD"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
